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Abstract 

Tidal dynamics of sliallow estuaries and lagoons is a complex matter that has attracted the attention of a large number of 
researchers over the last few decades. The main purpose of the present work is to study the intricate tidal dynamics of the 
Tagus estuary, which states as the largest estuary of the Iberian Peninsula and one of the most important wetlands in 
Portugal and Europe. Tagus has large areas of low depth and a remarkable geomorphology, both determining the complex 
propagation of tidal waves along the estuary of unknown manner. A non-linear two-dimensional vertically integrated 
hydrodynamic model was considered to be adequate to simulate its hydrodynamics and an application developed from the 
SIIV1SYS2D model was applied to study the tidal propagation along the estuary. The implementation and calibration of this 
model revealed its accuracy to predict tidal properties along the entire system. Several model runs enabled the analysis of 
the local variations in tidal dynamics, through the interpretation of amplitude and phase patterns of the main tidal 
constituents, tidal asymmetry, tidal ellipses, form factor and tidal dissipation. Results show that Tagus estuary tidal dynamics 
is extremely dependent on an estuarine resonance mode for the semi-diurnal constituents that induce important tidal 
characteristics. Besides, the estuarine coastline features and topography determines the changes in tidal propagation along 
the estuary, which therefore result essentially from a balance between convergence/divergence and friction and advection 
effects, besides the resonance effects. 
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Introduction 

The water movements and the turbulent mixture that result 
from the tidal forcing express problems and interesting challenges 
in the hydrodynamic field. Knowledge of tidal heights and tidal 
currents structure is therefore essential to understand problems 
such as dispersion rate of pollutants, sediment transport and 
erosion processes in coastal areas [1]. Moreover, tidal asymme- 
tries strongly influence nutrient balances, sediment loads, 
particles and pollutants transportations, etc. [2]. Thus, the 
understanding of the central processes lined by the tidal wave 
propagation seems to be crucial to obtain an overview 
concerning to the different uses of the coastal systems and in 
this particular study, of the Tagus estuary (Fig. 1), which 
dynamics is mainly ruled by the tidal forcing [3]. 

This coastal system is one of the largest estuaries of the West 
coast of Europe and is located in the most populated area of 
Portugal, including the capital Lisboa [4] . Likewise, Tagus estuary 
represents an important ecological system that is classified as a 
natural reserve since 1976 (Ramsar Convention on Wetlands). 
However, most of its margins are occupied by anthropogenic 
activities [5], and consequently is stUl subject to considerable 
urban and industrial pressures that induce several environmental 
problems including contamination by heavy metals and faecal 
material. In fact, this estuary plays an important role at a 



metropolitan and national level, due to two main aspects: it has a 
central environmental element in the Lisboa Metropolitan Area, 
which arises from its biological potential and its natural and 
vmique heritage richness; it concentrates several uses and 
economic and social functions - tourism, recreation, leisure, 
industry, agriculture, aquaculture, ports, urban areas - which need 
to be in balance with development [6]. Nevertheless, both the 
ecological and human values may be affected by near future sea 
level rise (SLR) [5]. 

[7], for example, in a study about SLR impacts in Tagus estuary 
salt marshes dynamics, indicates that instantaneous velocity 
intensity will increase in a climate change scenario, promoting a 
rising tendency for the sediments exportation outwards the system, 
which represent a serious threat for the stability of the system 
margins and consequentiy for the surrounding ecosystems. These 
expected hydrodynamic changes wUl increase margins instability 
and contribute to coastal erosion. Furthermore [8], in a study 
about the vulnerability of the Tagus estuary margins, indicate an 
increasing vulnerability of urban areas to SLR, being important 
the knowledge of the present system patterns to develop future 
planning, management guidelines and adaptive measures to SLR. 
Deeply studies about SLR effect in Tagus estuary hydrodynamics 
will be crucial to develop strategies to minimize climate change 
impacts. To support such studies, it is essential to understand the 
current hydrodynamic features and patterns of the Tagus estuary. 
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Figure 1. Geographic location of Tagus estuary and the numerical bathymetry, with depth in meters relative to the chart datum 
(2.0 m below the mean sea level] and with locations of the stations used in the numerical model calibration. 

doi:1 0.1 371/journal.pone.0080450.g001 



However, numerous studies have been performed concerning 
the biogeochemical dynamics of this estuary ([9], [10], [11], [12], 
[13] and [14]) while despite its relevance in the overall estuarine 
dynamics, the amount of work dedicated to its physical dynamics, 
and especially to its tidal properties, is limited. The contributions 
found in this field are in fact restricted. [15] and [16] performed a 
numerical study considering the 3-dimensional currents dynamics 
at the estuary mouth and the effect of tidal flats on estuary 
hydrodynamics. Recently [3], described the main forcing mech- 
anisms role on the estuarine dynamics, finding that tide is the 
estuary key forcing mechanism, although river discharge influ- 
ences the water column stratification in a seasonal time scale. The 
meteorological effects (wind and air temperature) are essentially 
seasonal, influencing the surface water temperature [17]. investi- 
gated the Tagus estuary plume dispersion induced by wind and 
river flow forcing using a 3-dimensional nested models set, 
concluding that the plume dynamic besides the importance of 
the freshwater discharge is highly influenced by the coastline 
geography and by the wind regime [18]. performed a preliminary 
study of Tagus estuary tidal dynamics through an analysis of the 
amplitude, phase and tidal ellipses parameters of the main semi- 



diurnal and diurnal tidal constituents, revealing that they are 
extremely dependent on the estuarine topography and coastline 
features. 

This work aims to provide a deepest analysis of local variations 
of the Tagus estuary tidal dynamics, also giving support to ongoing 
biogeochemical research. More precisely, this work aims to deepen 
the previous study by [18], verifying the relation between the 
estuary dynamics and its topography and coastline characteristics, 
thus presenting further details and actual information about the 
Tagus estuary hydrodynamic features. The hydrodynamic model 
setup and calibration is also described in order to demonstrate the 
accuracy of this modelling application to study and understand the 
Tagus tidal dynamics. For that, a numerical modelling study of the 
tidal wave propagation along the estuary was performed, seeking 
to determine and to analyse the amplitude and phase patterns of 
the main tidal constituents, the tidal asymmetry, the tidal ellipses, 
the form factor and the tidal dissipation. Besides the Tagus 
dynamics knowledge improvement, these results will also contrib- 
ute for a deepest understanding of the biogeochemical processes in 
this system. 
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Study Area 

Covering an area of approximately 320 km^ and with a mean 
volume of 188 km^, the Tagus estuary is one of the largest 
estuaries of Europe and it is the most extensive wetland area of the 
Portuguese territory ([7] and [19]). Its width varies between 400 m 
at its head and 1 5 km in the central bay and has an average depth 
of 5,1 m [16]. 

Morphologically, Tagus can be divided in two distinct regions 
that present different morphologies and properties: the lower and 
upper estuary (Fig. 1). The lower estuary is a channel of about 
30 m depth that connects with the Atlantic Ocean and opens in a 
large bay (upper estuary) on the east side. The upper estuary (bay), 
widest and shallowest, extends from Vila Franca de Xira to the 
main channel (lower estuary) and it is characterized by extensive 
zones of tidal flats, salt marshes, small islands and a net of narrow 
channels [7] and [15]. 

The Tagus river is the major estuarine source of freshwater, 
with an annual average flow of ~370 m'^/s [3]. Other freshwater 
inputs, the Sorraia and the Trancao rivers, are comparatively 
small, with average annual discharges of ~35 and ~2.5 m^/s, 
respectively [3]. According to [3], the river discharge forcing 
mechanism acts essentially on a seasonal time scale, mainly 
influencing the surface layer salinity. In fact, Tagus has a mean 
tidal prism of ~60xl0'^ m'^ and the river inflow per tidal cycle is 
-8.2x10*" m'' [3], which means that the freshwater input is only 
about 1.4% of the marine water input during the same period. 
The estuary should be thus vertically homogeneous during large 
periods and dominated by tidal forcing. 

The Tagus is ebb dominated with floods typically one hour 
longer than ebbs (stronger \'el()cities during ebbs), and thus 
inducing a net export of sediments [16]. The ar(;a affected l)y tides 
reaches 80 km landward of Lisboa and maximum tidal currents 
achieve ~2.0 m/s [12]. This mesotidal estuary presents an 
average tidal range at the seaward end of 2.4 m, and an extreme 
range changing from 0.9 m in neap tide to 4.1 m in spring tide. 
Tides are primary semi-diurnal and the tide form number is equal 
to 0.10 in Cascais (below the limit for mixing tides) [15]. Af^ is the 
dominant tidal constituent with amplitudes of ~ 1 m. The phase 
difference between Cascais and Vila Franca de Xira is 
~1 h20 min for semi-diurnal constituents [16]. According to 
[16], the amplitudes of astronomic constituents grow rapidly in the 
lower estuary and more steadily in the upper estuary, up to Povoa 
Sta. Iria, and then decrease up to Vila Franca de Xira. The 
shallow water constituents are negligible in the coastal region, but 
increase rapidly inside the system [16]. 

Hydrodynamic Model 

Previous hydrodynamic modelling implementations for the 
Tagus estuary in general have considered barotropic models, 
since the estuary is well mixed or partially mixed, except in flood 
scenarios [3]. Therefore, in this work, a modified version of 
SIMSYS2D ([20], [21] and [22]), based on the shallow water 
equations is used, expressing the mass and momentum conserva- 
tion, and assuming incompressibility and hydrostatic pressure. In 
order to describe the relevant areas of tidal flats in this estuary, the 
model grid cells are allowed to dry/flood in a complex process 
desc:ribed in [20] . This model was recently used with success in the 
hydrodynamic study of several coastal areas, as: Ria de Aveiro 
lagoon (Portugal) [22], [23], [24], [25], [26] and [27]; Lima 
estuary (Portugal) [28] and [29]; Tagus estuary (Portugal) [18], 
Patos lagoon (Brazil) [30] and Maputo bay (Mozambique) [31]. 

In this particular case, before applying this model it is essential 
to describe the hydrodynamic model setup and caUbration for 



Tagus estuary, in order to demonstrate the accuracy of this 
modelling application to study and understand the Tagus tidal 
dynamics. 

Model setup 

The major limitation of previous modelling applications 
developed for Tagus estuary has been the low horizontal 
resolution of the grids used, insufficient to properly resolve all 
channels [32]. In this work we chose to develop a grid with a 
resolution higher than that used in previous studies, considering 
the dimensions Ax = A)i = 200 m. The numerical bathymetry 
(Fig. 1) was developed from topo-hydrographical data measured 
by the Hydrographical Institute of the Portuguese Navy. A 
Krigging method was used to generate the rectangular computa- 
tional grid, resulting in 227 cells in the x direction and 202 cells in 
the_>' direction. 

In order to optimize the computational time and to obtain the 
best results, assuring the model stability, several simulations were 
performed to determine the ideal time step and horizontal viscosity 
to use in the model. For fixed horizontal viscosity and time step, 
different values of time step and horizontal viscosity were adopted, 
respectively, and the model results for the stations represented in 
Figure 1 were analys("d. 

The best model results are achieved with time step and 
horizontal viscosity of 2 minutes and 15 m^s respectively. 
Several simulations were also performed to verify the model results 
independence relatively to the initial conditions, from which was 
concluded that the model spin-up period is about 48 hours. 
Consequentiy, these time step, horizontal viscosity and spin-up 
period values were included in all the study simulations. 

The sea surface elevation (SSE) at the ocean open boundary 
results of the harmonic synthesis values that were determined for 
Cascais, using harmonic constants for this tidal gauge. Once the 
present study intends to analyse the astronomic tide propagation 
and the open boundary is very close to the tidal gauge position, 
this procedure minimizes the model errors since it as the 
advantage of remove other signals besides tides from the observed 
signal such as wind and variations of atmospheric pressure. 

Model Calibration 

A mathematical model is by definition an attempt to 
approximate and reproduce real phenomena [33]. The approx- 
imations and parameterizations used for the model synthesis lead 
to discrepancies and deviations of model results from nature. 
Therefore, before being applied to a specific location, models 
should be calibrated and sometimes validated. Model calibration 
appears in various forms, depending on data availability, on water 
body characteristics, and most of all, on modellers' perceptions 
and opinions [34]. 

In this work, three distinct and complementary hydrodynamic 
model calibration methods are used, according to the methodol- 
ogy proposed by [22], and followed by [26], [35] and [36], in 
which an initial qualitative calibration is followed by the 
quantification of the models accuracy. The first step consists in a 
visual direct comparison between model results and observations, 
changing the model parameters in order to optimize their 
adjustment. In the second step, the root-mean srjuare (RMS) 
errors and the Skill assessment between the SSE model results and 
observed data are calculated. The third step compares the SSE 
harmonic constants of the observed data and model results. 
Harmonic constants for 1 3 tidal gauges are available for the model 
calibration, providing an excellent coverage of the estuary (Fig. 1). 

Comparison between predicted and observed series of 
SSE. Time series of SSE were synthesized for the 13 calibration 
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stations applying the t_tide matlab® harmonic analysis package 
[37] to the available observed data (which consists in harmonic 
constants). It was reproduced the period from March 28 to May 1 
of 2012. This period was also reproduced numerically, in order to 
obtain predicted time series for each one of the stations. Several 
simulations were performed, adjusting the Manning coefScient in 
order to optimize the adjustment between the model results and 
the observations. The parameter used in this adjustment during 
model calibration is the bottom stress represented in the Manning- 
Chezy formulation, being 0.04 the constant Manning coefficient 
value that induces model results that best fit the observations. 

Figure 2 shows the comparison between the predicted (lines) 
and the observed (points) time series of SSE, for 6 of the stations 
used in the model calibration. The final adjustment between the 
series is good both in amplitude and phase although, during 
spring tides, sea surface elevation is underestimated by the model 
in some stations. The fortnightly tidal cycle is also well 
represented in the predictions, with a good reproduction of 
amplitude and phase, both in neap and in spring tides. 

As the model open boundary condition is the SSE time series 
from Cascais, the agreement between the observed and predicted 
series for Cascais and Pa^o de Arcos is perfect, because these 
stations are very close to the ocean open boundary. The adoption 
of the observed SSE time series at Cascais as boundary condition 
allowed an optimization of the model results for the entire estuary, 
since obtaining a perfect match between predicted vs. observed 
time series at the lower estuary minimizes the errors for the 
remaining stations. 



RMS errors and Skill assessment. The next step in the 
model calibration is the quantification of the model accuracy 
through the determination of the deviations between the model 
predictions and the observations. With this purpose were 
computed the root-mean square (RMS) errors and the Skill. 

The RMS error of the difference between the predicted and 
observed SSE series was previously used by several authors (for 
example [15], [22] and [25]) to evaluate model's accuracy; 

fi y 

^^^=i]^E^"o(?,-)-c.w] I (1) 

where Coft) and Cpft) are the observed and the predicted SSE, 
respectively, and Mo is the number of measurements in the time 
series. According to [38], the RMS values should be compared 
with the local tidal range and if they are less than 5 % of the local 
amplitude, the correspondence between the model results and 
observations is considered excellent. If these values vary between 
5% and 10% of the local amplitude, the match is considered very 
good. 

The model predictive Sldll is being recentiy used to evaluate all 
prognostic quantities and is based on the quantitative agreement 
between model results and observations. This statistical method 
was developed by [39] and recently used by [38], [40] and [41]. 
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Figure 2. Comparison between the predicted (lines) and the observed (points) sea surface elevation time series for the stations of 
Cascais, Pa^o de Arcos, Trafaria, Cacilhas, Lisboa, and Montijo. 

doi:10.1371/journal.pone.0080450.g002 
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Skill = 1 _ , (2) 

E(IC/,-Col + ICo-Col) 

A Skill value equal to one corresponds to a perfect agreement 
between model results and observations while a complete 
disagreement yields a Skill of zero. 

The RMS and Skill were computed for each calibration station 
and are presented in Table 1. The analyses of these results reveals 
that the model accurately reproduces the tidal wave propagation 
along Tagus estuary, since all the Skill values are very close to 1 
and, likewise, RMS values are also between 0% and 15% of the 
local amplitude for almost stations. 

The poorer agreement was found for Seixal, Povoa de Sta Iria, 
Ponta da Erva and Alcochete stations, because RMS and Skill 
values do not account the diflference between the predicted and 
observed series in low tide. These stations are located in salt marsh 
zones, where the model simulates the wetting/ drying along typical 
tidal cycle of tidal flats, with null SSE at the low tide. As the 
observed series were determined through harmonic synthesis, this 
effect is not observed in the data for these stations, and SSE values 
different from zero are predicted for low tide. 

Harmonic analysis. Comparative harmonic analysis of 
model predictions and observations is another method widely 
used to evaluate the accuracy of the model's performance ([22], 
[25], [26] and [35]). This method has the advantage of allowing 
the quantification of the models performance for each of the main 
tidal constituents, according to the tidal wave propagation 
characteristics along the estuary. In this study the harmonic 
analysis was performed using t_tide matlab® package [37] 
considering the harmonic constants values that were calculated 
from the model predictions for the 13 stations, which allowed the 
determination of the amplitude and phase of the main constitu- 
ents, as well as the respective errors (Fig. 3). The period in analysis 
is the same that was used in the previous comparisons. Only the 
main harmonic constituents are presented, since the entire 
harmonic constituents list is extremely extensive. The predicted 
amplitudes interval for the semi-diurnal constituents (A/^, and 



Table 1. RMS and predictive skill for al 
stations. 


the calibration 




station 


RMS (m) 


SKILL 


Cascais 


0.015 


1.000 


Pa^o de Arcos 


0.043 


1.000 


Trafaria 


0.117 


0.999 


Cacilhas 


0.119 


0.999 


Lisboa 


0.148 


0.999 


Seixal 


0.634 


0.980 


Alfeite 


0.289 


0.996 


Montijo 


0.255 


0.997 


Cabo do Ruivo 


0.231 


0.997 


Alcochete 


0.377 


0.993 


Ponta da Erva 


0.705 


0.974 


Povoa de Sta. Iria 


0.358 


0.994 


Vila Franca de Xira 


0.234 


0.997 


doi:1 0.1 371 /journal.pone.0080450.t001 



jV^) includes the observed amplitudes at most stations, being 
the only exception for stations farther away from the estuary 
mouth (Ponta da Erva, Povoa de Sta Iria and Vila Franca de 
Xira). For the phase, the results were good in all the stations, 
except a small deviation found for jV^ in Trafaria. For the diurnal 
constituents {Kj and 0/), the observed amplitudes met inside of the 
predicted interval for all the stations. For the diurnal phase the 
results were also very good since the only exceptions were found in 
Trafaria, Alfeite and VUa Franca de Xira (for Kj) and in Povoa de 
Sta Iria (for Oj). For the lunisolar synodic fortnightly constituent 
MS( the results were good in all stations both in phase and 
amplitude. For the shallow waters constituents (M4 and M(j the 
predicted interval includes observed amplitudes in most stations, 
except in Lisboa, Montijo and Vila Franca de Xira for M4. 
Concerning the phase, the results were good only for Cascais and 
Pago de Arcos stations for M4 and Cascais and Montijo for Mg. In 
the remaining stations there was a great difference between 
observed and predicted phases, but it should be enhanced the 
great errors associated with the determination of the predicted 
phases for shallow water constituents. 

The significant errors found for the shallow water constituents 
are not unexpected, as overtides are usually more difficult to 
predict correctly than the major astronomic constituents. M4 and 
Affi prediction errors may be due to bathymetric inaccuracies that 
cannot be corrected by the adjustment of the bottom friction 
coefficient. Even if these discrepancies are partially mitigated by 
the calibration for the major constituents, friction influences 
differently the non-linear constituents. However, as these constit- 
uents are generated through the advective and finite amplitude 
terms of the shallow water equations, which are found to be 
accurately predicted by the model, the results are expected to be 
precise enough to reproduce tidal asymmetry in this study. 

Summary. The results show that the amplitude and phase of 
the tidal constituents can be considered satisfactorily represented 
by the model for the entire estuary. The deviations found between 
model predictions and observed data are not significant for most 
stations and constituents. Considering the direct comparison of the 
observed and predicted series of SSE, as well as the RMS and SkUl 
values, it is verified that the best model results are found close to 
the estuary mouth but, generally, there is also a good adjustment 
for all the stations situated both in the central zone and far from 
the estuary mouth. The results were less satisfactory for stations 
close to the estuary margins, where the bathymetry is further 
complex (narrow channels and small bays) and the water is shallow 
(salt marsh and tidal flat areas). These higher discrepancies may be 
due to the inexact definition of the bathymetry (due its complexity, 
it is difficult to numerically define it) or to the natural difficulty in 
reproducing accurately the wetting/ drying of these shallow areas. 

It should be pointed out that the skill and RMS values are 
similar to those determined by [42] for the Guadiana estuary, [30] 
and [43] for the Patos lagoon, [40] for the Hudson River estuary, 
[22], [25], [26] and [35] for the Ria de Aveiro lagoon or [38] for 
Ria Formosa lagoon. The hydrodynamics of these estuarine 
environments was studied using finite elements and finite 
diflerence models. Those authors have found predictive skills 
ranging from 0.90 to 1 and RMS values of tidal elevations around 
10% of the local amplitudes. As such, the results obtained show 
that the 2D application developed for Tagus estuary reproduces 
accurately the propagation of the tidal wave and thus the 
barotropic flows in this complex estuary. Therefore, it was 
concluded that this 2D model can be considered successfully 
calibrated and thus can be used to study important issues 
concerning the local variations in the estuary tidal regime. 
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Figure 3. Comparison between the main harmonic constants amplitudes (a) and phases (b) determined from predicted (red) and 
the observed (blacl<) sea surface elevation time series for the stations of Cascais, Pa^o de Arcos, Trafaria, Cacilhas, Lisboa, and 
Montijo. The values from predictions include error bars. 
doi:1 0.1 371 /journal.pone.0080450.g003 



Tidal Characterization of Tagus Estuary 

Due to the small dimensions and shallowness of Tagus estuary, 
its tidal properties are essentially driven by the oceanic tidal 
characteristics at its entrance and by changes and co-oscillation 
that results from tidal propagation along the estuary. Consequent- 
ly, the analysis of changes in astronomical constituents along the 
estuary becomes essential to understand the factors ruling its tidal 
dynamics. 



In the present work, the astronomical constituents and Oj 
wiU be analysed because previous studies revealed that they are the 
major semidiurnal and diurnal constituents, which was also 
verified in the previous section results. M4 is the most important 
shallow water constituent and has relevance in the definition of the 
local tidal asymmetry, and A/g is important in the bottom friction 
effects analysis, and consequently will be also analysed. 

The present section includes the results of the numerical 
simulations that were carried out in order to establish the tidal 
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characteristics of Tagus estuary. The parameters defined during 
calibration of the model were kept and a time series of sea surface 
elevation (SSE) was applied for the period from October 8 (201 1) 
to November 13 (2012) as ocean boundar)' condition. Hourly 
temporal estimation of SSE and horizontal fields of SEE 
throughout the computational domain were obtained. The time 
series of the velocity components were also computed at all grid 
points. The phase and amplitude distribution of the main 
harmonic Afj, 0;, M,y and M(; were determined and analysed 
based on those results by applying t_tide matlab® subroutines [37]. 
From this data set the form factor (requiring also the determina- 
tion of amplitudes distribution for and Ki), tidal ellipse, tidal 
asymmetry and tidal dissipation were also determined. 

Amplitudes and phase patterns of the main tidal 
constituents 

As general pattern, it was observed that the variations in 
amplitude occur where changes in the estuary morphology and 
depth are larger, what happens in the transition zone between the 
main channel and in the initial and upper bay area. The results 
also reveal that the amplitude of the tidal diurnal constituent is one 
order of magnitude lower than the semi-diurnal in the whole 
system. 

Concerning the amplitude of the higher harmonic constituent, 
M2, it increases progressively along the estuary until North of 
Cabo Ruivo (middle bay), diminishes until Povoa de Sta. Iria, and 
raise again until Vila Franca de Xira (upper estuar)'). It ranges 
from 0.97 m at the mouth of the estuary to about 1.11 m near Vila 
Franca de Xira. The tidal amplification found at the middle bay is 
unexpected, since usually occurs a decrease in the amplitudes of 
the principal diurnal and semidiurnal constituents as tidal wave 
propagates along the most estuarine systems, due to friction and to 
non-linear transfers of energy [44]. However, as the tidal wave 
propagates into shallower waters within the estuary, it becomes 
distorted due to several physical processes, including reflection at 
the estuary head, which is responsible for standing-wave gener- 
ation and local resonances [45]. Thus, the amplification found 
may be explained considering an estuarine resonance mode with a 
period close to twelve hours, due to the estuary length, which 
enhances the M2 constituent as suggested by [16]. Following [46] 
and [47] quarter-wave resonance in a tidal basin occurs when: 



where L is the Ic-ngtli of the basin, w is the tidal period and c is the 
phase speed of the tide. For the cxmstituent, 
a)=1.41xl0 *s 'and the mean water depth determined from 
the numerical bathymetry developed for Tagus estuary in the 
frame of this study is 5.1 m, resulting in the phase speed, 
c=7.07ms~'. Thus, for quarter-wave tidal resonance, 
Z.= 78.76 km which is very close to the actual 80 km length of 
Tagus estuary. This result is therefore consistent with the tidal 
amplification found along the estuary for M2 constituent. 

Other authors have also reported this local resonance. 
According to [48], tidal amplitude inside the Tagus estuary is 
larger than at its mouth as a result of a small resonance effect and 
therefore he considered that this estuary can be hydrodynamically 
interpreted as a greater tidal basin connected to the ocean by a 
relatively straight and narrow channel. [15], for example, referred 
that tidal form number decreases along the system due to a 
resonant mode that increases the amplitudes of semidiurnal waves 
by roughly 40% in the upper estuary, leaving the diurnal waves 
mostiy unchanged. 



The tidal amplitude decline from Cabo Ruivo to Povoa de Sta. 
Iria is due to the shallow depths of the upper bay that induce 
higher dissipation of tidal energy- by bottom stress, therefore 
reducing the tidal amplitude. Upstream Povoa de Sta. Iria the 
channel convergence increases again the tidal amplitude. In some 
narrow and shallow channels as Serxal and Ponta de Erva M2 
amplitude decreases to 0.50 m as result of the large bottom friction 
influence. 

Regarding to phase, difierences in tidal propagation were also 
found along the two main areas of the estuary, the upper and 
lowc-r estuary, with M2 velocity being strongly attenuated in the 
shallow areas. In fact, it was observed a minor phase delay as tidal 
wave progresses toward the inshore zone of the system along the 
deepest lower estuary. It was also observed that the tide 
propagation delay is further significant in the bay due to the 
estuary morphology divergence and to depth decrease, which, in 
turn, makes bottom friction plus relevant inducing a lower tidal 
propagation velocity. Figure 4 shows that M2 phase is ~73° near 
the estuary mouth and ~125° in Vila Franca de Xira, which 
means that the propagation delay for this constituent is about 
1 hour and 47 minutes (52°). This value is slightly higher than that 
presented by [16], which refers to a delay of 1 h and 20 minutes. 

The diurnal constituent 0/ has amplitude of ~0.060 m near to 
the mouth, which keeps relatively constant along the estuary, 
reaching ~0.064 m near to Vila Franca de Xira. Once again, in 
the channels mentioned above (Seixal, Ponta de Erva and Povoa 
de Sta. Iria), the amplitude of the constituents drops to ~0.050 m. 

With respect to phase lag, Oi shows a delay of about 1 hour and 
50 minutes, which corresponds to a phase difference of 25°: near 
the mouth of the estuary is close to 325° and at the head is about 
350° (Fig. 4). 

Amplitudes and phase pattems of the shallow water 
constituents 

The sea surface water elevation (Q is usually a non-negligible 
fraction of the total water depth in estuarine shallow waters. As 
such, it has to be considered in the tidal wave celerity: 

c = {pg{HH)f'^ (4) 

where p is water density, g the gravity acceleration and H the 
water depth. The wave crest travels rapidly than the troughs, 
resulting in shorter and further intense flood than ebb currents. 
Bottom friction also influences tidal wave propagation in shallow 
waters, delaying the water movement. As these processes depend 
on the- square of current velocity, overtide generation is non-linear. 
The shallow water M4 or Mg constituents, both originated from 
M2, are therefore overtides that represent these distortions of the 
tidal levels normal harmonic variations, induced by shallow-water 
and friction effects. M4 is essentially generated by the non-linear 
terms in the equations of continuity and momentum (advection), 
due to different celerity of crests and troughs, but friction can also 
play a minor role [49] . Mg is basically generated by the quadratic 
friction term. 

The results for shallow water constituents demonstrate that their 
amplitudes are negligible in the coastal region, but increase rapidly 
inside the estuary (Fig. 5), being M4 the most significant of these 
constituents. As a shallow water constituent, M4 is strongly 
influenced by changes in system morphology and bathymetry. 
Therefore, amplitude and phase variations occur in areas where 
significant changes on system morpholog)- and bathymetry arise. 
Figure 5 shows that M4 systematically grows from Cascais to Cabo 
Ruivo, being the gradient higher between CaciUias and Cabo 
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Figure 4. Semi-diurnal M2 and diurnal constituents amplitude and phase for Tagus estuary. 

doi:10.1371/journal.pone.0080450.g004 



Ruivo. From this point on, its amplitude decreases upstream 
towards Povoa de Sta. Iria, and increases again with a strong 
gradient until Vila Franca de Xira. The amplitude varies from 
0.02 near to the estuary mouth to 0.05 in the upper bay area, and 
close to Vila Franca de Xira the M4 amplitude is ~0.1 m. In the 
channel nearby Montijo it is observed a significant increase in 
amplitude from 0.06 m to 0.08, which can be explained by the 
shortening and shallowness features of this area. It is important to 
point out that in the shallow inner bay (between Lisboa and Vila 
Franca de Xira), M4 has higher amplitude than the diurnal has. 

The Mg amplitude increases upstream, weakly until Cabo Ruivo 
and then with a stronger gradient toward Povoa de Sta. Iria, where 
it reaches its maximum values. In the main channel, the amplitude 
is relatively constant, approximately 0.013 m. However, from the 
initial area of the bay to the upper bay, is observed an amplitude 
increase of ~0.018 m: the amplitude varies from ~0.014 to 
~0.032 m, respectively. Near Vila Franca de Xira, Figure 5 shows 
a decrease of ~0.005 m, which means that in this area the 
amplitude value is about ~0.027 m. 

As Cabo Ruivo and Povoa de Sta. Iria seem to be turning points 
for the shallow-water constituents, this might correspond to 
changes in the respective generating processes. [16] also reported 
this fact, stressing the dominance of advective accelerations 
(responsible for M4 generation) in the deeper part of the estuary 



compared with the larger influence of friction (responsible for A'lg 
generation) in the shallower part of the estuary (Fig. 1). 

Regarding to the M4 phase lag, in the initial zone of the estuary, 
the propagation delay is very small: from the mouth to the upper 
bay area, the phase difference is about 15°. However, this delay 
increases upstream and the difference becomes more significant 
near Vila Franca de Xira (varies from ~60° to ~210°, which 
corresponds to 2 hours and 30 minutes) as a result of strong 
irregularities in the morphology. Concerning to Mg phase, the 
lower values (~75°) are observed at the upper bay, characterized 
by the shallower depths, which induce the higher bottom friction 
and where the Mg higher amplitudes were found. As a result it 
should be concluded that this constituent is essentially generated in 
this area, propagating up and downstream along the estuary. The 
phase increases to ~ 1 1 4° at the end of the main channel and to 
~160° at the mouth, and to ~124° near Vila Franca de Xira. 

Form factor 

The relative importance of diurnal and semi-diurnal constitu- 
ents could be expressed in terms of a form factor or number factor. 
In this study, the form factor for the entire estuary was assessed, 
derived from the harmonic analysis results and applying the 
following expression [50]: 
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Figure 5. Shallow-water M4 and Me constituents amplitude and phase for Tagus estuary. 

doi:10.1371/journal.pone.0080450.g005 



M2 + S2 



(5) 



where Oy, Kj, and S2 are the amplitudes in meters of 
correspondent constituents. Tide is semidiurnal if 0<i^<0.25; 
mixed tide - semidiurnal tides predominate if 0.25<F<1.5; mixed 
tide - diurnal tides predominate if 1.5<F<3.0; and, finally, the 
tide is diurnal if ^>3.0 [45]. 

As this estuary is mesotidal and tides are primary semi-diurnal, 
M2 and S2 amplitudes are significantly higher than those of diurnal 
constituents Oj and Kj, being thus expected a low form factor F. In 
fact, Figure 6 demonstrates that F value is less than 0.25 
throughout the estuary, meaning that tide is semidiurnal. In fact, 
F is approximately constant and equal to 0. 1 within the system. 
Although almost imperceptive due to the colour scale that was 
necessary to applied to represent the model results, it is verified a 
slightly decrease of the form factor along the system, mainly within 
bay, mostly as a result of M2 increased amplitude, as earlier 
discussed. Tidal form number decrement along the system was 
described by other authors, such as [15], referring that this pattern 
results from the resonant mode that increases the semi-diurnal 



waves amplitudes inside the estuary. In shallowest areas such as in 
the upper bay margins and mostly in channels along the Sebcal 
and Povoa de Sta. Iria, form factor increases, approaching 0.49, 
due to the strong decrease of the semi-diurnal constituents 
amplitudes in these areas, as was previous noted. 

Tidal Ellipses 

The determination of tide ellipses for the main tidal constituents 
allows the assessment of tidal currents patterns along the estuary 
[31]. Figure 7 shows the tidal ellipses for semi-diurnal, diurnal and 
shallow water constituents in the Tagus estuary, and from its 
analysis is found that in general all ellipses have high eccentricity, 
meaning that tidal currents are mainly rectilinear and the ellipses 
sense of rotation should not be significant. This pattern is induced 
by the constraints imposed by local bottom topography and 
geometry. Therefore, the ellipses orientation (for semidiurnal, 
diurnal and Mg constituents is similar) reveals the importance of 
these factors on tidal fluxes, with main flux occurring essentially 
through deeper zones. In fact, tidal ellipses results revealed 
preferential flow channels in the bay, with highest speeds being 
reached in two deepest areas. From Figure 7 it is also observed 
that in the Tagus estuary currents attain maximum values at the 
system mouth (~0.8 for the M2 constituent; ~0.035 for Oj 
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Figure 6. Form factor distribution in the Tagus estuary. 

doi:1 0.1 371 /journai.pone.0080450.g006 

constituent). The general patterns found for both constituents 
show that current velocity decreases from the lower to the upper 
estuary. Lower current values were found in the bay, near to land 
boundaries. In the upper zone of the bay, the current loses 
intensity and values decrease, approaching 0.3 m/s in the case of 
M2 and approximately 0.02 m/s for 0/. 

In turn, M4 shows a difiFerent pattern. As a shallow water 
constituent, the highest values are found when significant 
differences in morphology and bathymetry arise (upper zone of 
the bay, essentially), being ~0.08 m/s. In the main channel, 
velocity is about 0.04 m/ s and in the initial zone of the bay its 
values decrease to 0.02 m/s. It should be referred that M4 current 
intensity is locally important, since it is often larger that 0; 
intensity, making this contribution the 4* most important in 
Tagus tidal dynamics. The pattern found for Mg is similar to those 
found for semi-diurnal and diurnal constituents, although with 
very low current intensities. 

Tidal Asymmetry 

A symmetrical tide occurs when the tide rise and fall have 
identical duration, presenting approximately equal maximum 
velocities, which results in no net overall sediment transport. Tidal 
asymmetry happens when the ebb and flood durations are unequal 
and it is caused by tidal wave distortion during propagation into 
shallow waters, along the coastal shelf and on entry into estuaries. 
This phenomenon is induced by nonlinear effect mechanism of 
tidal propagation, resulting from finite amplitude effects on the 
non-linear advection and friction terms. Likewise, information of 
tidal asymmetry may become very useful to determine currents 
patterns, including the current flood or ebb prevalence [51]. In 
many coastal worldwide areas where tide is semi-diurnal, as in the 
Tagus estuary, the dominant tidal constituent is and thus the 
most significant shallow water constituent is M4, the first harmonic 
of Afj. However, as observed in section 4.2, Mg should also be 
considered in the case of Tagus estuary. Tidal distortion 
(asymmetry) can be represented by harmonics of astronomical 
tidal constituents [44] . Thus, the ratios of M4 and and of Mg 
and M2 amplitude (which corresponds to the A, values) can be 
analysed to determinate the tidal asymmetry magnitude generated 
within the Tagus estuary. Similarly, relative phases of M2 and M4 
and of M2 and Mg determine the type of asymmetry. Therefore, 
the asymmetrical coefficients (amplitude ratio) and (p (relative 



phase) can be represented by the following equations: 

M2 M2 



cp =26,. -6,. ;(p,. =26,. -6,.. (7) 

where A is the amplitude, 6 indicates the phase and Mg, M4 and 
M2 correspond to the tidal constituents. The flow is flood dominant if 
0°<(p<180° and ebh dominant if 180°<cp<360°. Eventually, this 
creates net sediment transport over a tidal cycle, with flood 
dominated currents resulting in net sediment transport into the 
estuary and, in the other side, under ehb dominant conditions, the 
opposite happens and net seaward transport will take place 
causing sediment export from the estuary. In fact, tidal asymmetry 
is frequently the dominant factor in causing net sediment transport 
and deposition which results in sediment trapping in coastal areas 
and estuaries [52]. The navigability of estuarine channels and 
estuaries geological evolution are consequently affected by tidal 
asymmetry [44] which is, therefore, a fvmdamental factor for 
morphological development of tidal basins. 

Amplitude ratio and relative phase difference: M4. The 
analysis of Figure 8 reveals that there is a progressive increase of^^ 
along the estuary. Entering in the system, tide suffers a severe 
distortion as a result of the strong shortening that is felt in the input 
channel (lower estuary): from the ocean to the initial zone of the 
channel, amplitude ratio increases from 0.025 to 0.041, approx- 
imately. At the bay, where depth is less. A,- value exceeds 0.048. 
The dark colours areas (dark brown) northeast from the bay 
correspond to narrow and shallow channels and tidal flats, which 
promote energy dissipation due to the bottom friction effects and, 
consequentiy, represent areas with strong tidal distortion. Some 
deeper channels that do not uncover with the tide, as the Seixal 
and Povoa de Sta. Iria channels, present maximal asymmetry 
values, close to 0.549, as a result of the water velocity increment 
due to the channel convergence. 

Throughout the system, the amplitude ratio A^ is always higher 
than 0.02 (excluding the mouth of the estuary), which allows to 
state that tidal asymmetry occurs in the whole system. An 
asymmetric Tagus estuary behaviour is also mentioned by [3]. 

In some other coastal systems, such as the Ria Formosa [53], 
there are areas where amplitude ratio value is zero, so it does not 
occur tidal asymmetry, and others that reach values of 0.9. Given 
this example, is verified that tidal asymmetry in the Tagus estuary 
is relatively low, with values of A,- higher than 0.9 only in the 
narrower and shallower channels (as near Serxal and Vila Franca 
de Xira). 

For relative phase is concluded that results are in agreement 
with bibliography that state that this estuary is dominated by the 
ebb, which means that rising tides duration exceeds those of falling 
tides, causing a tendency for stronger ebb than flood tidal currents. 
Actually, Figure 8 shows that for the entire estuary (excluding 
narrower and shallower channels) relative phase value is over than 
180°, indicating ebb dominance. Relative phase values are maximum 
in mid-estuarine region (transition between the inlet channel and 
the inner bay) revealing that the highest distortion on water level 
variation towards ebb dominance occurs in this area. The ebb 
dominance has a strong impact on sediment dynamics, promoting 
sediment transport out of the estuary and therefore inducing 
erosion in these areas. However, in narrower and shallower 
channels relative phase differs from the rest of the estuary, 
similarity to the pattern observed for the amplitude ratio: relative 
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Figure 7. Tidal ellipses distribution of tidal constituents M^, 0„ M4 and Mg in Tagus estuary. 

doi:1 0.1 371 /journal.pone.0080450.g007 



phase is less than 180° (Seixal, Ponta de Erva and Vila Franca de 
Xira channels - Figure 9). Thus, Jlood dominance occurs, favouring 
sedimentation processes in these places with consequent accretion. 
Amplitude ratio and relative phase difference: Mg. The 

majority of studies on tidal asymmetry only analyse the 
importance of M4 on the water level distortion, thus only a 
small attention has been devoted to research contribution to 
tidal asymmetry. In this work, amplitude ratios AMg/AM^ and 
relative phase ^OM^-OMg were also determined from the 
harmonic analysis results, in order to identify a possible relation 
between these parameters and tidal asymmetry. Results for water 
level time series are presented in Figure 8, showing a progressive 
increase of A,- along the estuary. Although A, values are 
approximately constant in the main channel (about 0.013), at 
the bay (where the depth is lower comparatively to the main 
channel) values are significantly higher. In fact, there is an 
increase of 0.019 within the bay: A^ rises from 0.013 to 0.033 
from the initial to the upper bay. The large increase of A,- along 
the bay reveals the intensification of the relative importance of 
bottom friction, resulting in a low water delay [49] , and therefore 
contributing to the decline of the time between low and high 
water. Thus, A^ pattern at the upper bay can be probably related 
to flood dominance trend, previously identified though relative phase 



26M2-6M4. A, amplitude ratio is always higher than 0.01 
throughout the system, revealing that tidal asymmetry occurs in 
the whole estuary. 

For relative phase, excluding the initial zone of the main 
estuary where relative phase values of 152° are observed, values 
are always higher than 180° in the whole estuary. At the end of 
the main channel the relative phase values are about 190° and in 
the upper bay reach 342°, which means a difference of 152° 
along the bay. Close to Vila Franca de Xira relative phase values 
decrease to 338°. 

Tidal Dissipation 

In this study the mean rate dissipation of tidal energy per unit 
area due to bottom friction is estimated from the application of the 
following equation [47]: 



£= — 
T 



(8) 



where Co is the bottom friction used in the calibration of the 
model and U and Fare the x andjc components of the current, 
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Figure 8. Amplitude ratio {A^ and relative phase (Cp) for the and constituents in Tagus estuary. 

doi:1 0.1 371 /journal.pone.0080450.g008 
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Figure 9. Spatial-logarithmic distribution of the average dissipation of tidal energy in neap and spring tides in Tagus estuary. 
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respectively. The integration presented in equation (8) was taken 
over one tidal cycle for neap and spring tide cases in this study. 

Figure 9 shows spatial-logarithmic distribution of average tidal 
energy dissipation in neap and spring tides in the Tagus estuary. 
Energy dissipation values were analysed through e logarithm for a 
better interpretation of r("sults. 

According to [7], values of the average rate of energy 
dissipation arise, in general, associated with areas where 
boundary friction is high or where a sudden change of estuary 
geometry occurs. Thus, changes in estuary bathymetry and 
geometry could have a substantial effect in the way that tidal 
wave energy is dissipated and may even be responsible for tidal 
wave distortion. From Figure 9 is observed that the highest 
values are reached in the transition from the ocean to the input 
channel, in the main channel (lower estuary) and in the central 
deepest area of the bay (upper estuary) (Fig. 1). During neap 
tides, the value of maximal dissipation arises near the estuary 
mouth, in the narrowest part of the main channel and in the 
deepest area of the bay where depth varies between 0 and 12 
meters. Without considering the wetting/drying areas, the 
remaining places in the estuary- where maximal dissipation is 
found show that dissipation decreases with decreasing depth: at 
higher depths dissipation is close to 0.5 W/m^ and decreases 
progressively with depth, reaching values of ~0.2 W/m^. 

During spring tides, in the transition zone, in the main channel 
and in the central area of the bay, the maximal dissipation is 
~1 W/m^. Indeed, at these locations, dissipation is high and 
almost constant, except in the broad area of the channel where, 
along the shores, the dissipation is lower. In the bay, the energy 
dissipation decreases to ~0.03 W/m^ in shallower channels 
(depths between 0 and 2 meters). The closer shore areas 
corresponding to wetting/ drying zones have lower values. In the 
entrance channel dissipation is maximal, which can be justified by 
the considerable depths in the central area (turning the bottom 
friction negligible) and by the channel convergence and shortening 
(favours an increase in current velocity). Likewise, as the channel is 
narrow and the depth decreases along its banks, the boundary 
friction can also foster an increase in energy dissipation. In the 
bay, there is a strong depth decline comparatively to the entrance 
channel and a significant decline as well from the central area to 
the peripheral areas, resulting in a speed decrease in these areas. 
As in main channel, the central area of the bay is also a place 
where currents are strongest relatively to the peripheral areas that 
have shallower depths, which suggests that dissipation is associated 
with high current velocities. In areas of minor depth a decreased 
flow occurs and dissipation is actually smaller. Moreover, since the 
dissipation is apparently related to high current speeds, looking at 
Figure 9 and taking into account the maximum dissipation, it can 
be also concluded that flows occurs through preferential channels 
(as mentioned above), where depths are greater. 

The results show that tidal dissipation is higher during spring 
tides once tidal cunx-nts are weaker in neap tides, inducing less 
energy dissipation. In general, the maximum values of dissipation 
for neap tides are close to 0.5 W/m^ while for spring tides they are 
~ 1 W/ m^, showing that maximal dissipation during neap tides is 
about half of the maximum that occurs in spring tides. This 
pattern is in accordance with the results obtained from other 
coastal systems such as the Ria Formosa [53], Shark Bay [47], 
Maputo Bay [54] and Ria de Aveiro [55], where the tidal 
dissipation calculated for neap tides is considerable lower than for 
spring tides. Comparatively to Ria Formosa and Shark Bay, the 
results suggest that in general Tagus estuary represents a similar 
potential of energy dissipation once the maximum values in those 
coastal systems are approximately 1 and 0.7 W/m^. Maputo bay 



presents values of ~0.1 W/m over a large proportion of the bay 
although it also could reaches maximum values of 1 W/ m^ while 
it was identified maximum values of energy dissipation of 

~ 1 . 7 W/ m^ in Ria de Aveiro lagoon. 

The results of the present study are in agreement with that 
observed for the other studied parameters: Tagus estuary 
bathymetry, coasdine features and bottom friction have an 
extremely important role in the circulation patterns of this system. 
Previously authors refer this fact and, for example, [3] says that the 
Tagus estuary bottom topography is the main constraint to the 
tidal currents direction. 

Conclusions 

A two dimensional numerical model has been used to examine 
the tidal dynamics of Tagus estuary. The first step of this study 
revealed that the model accurately reproduces the tidal wave 
propagation in this estuary, so it can be considered successfully 
calibrated. Therefore, the tidal dynamics of Tagus estuary has 
been successfully resolved through the application of the 
SIMSYS2D model. It should be pointed out that the use of a 
2D model assures a good reproduction of the estuarine tidal 
dynamics, even neglecting the 3D baroclinic effects. Model results 
indicate that estuary- bathymetry and coastline features, namely 
the differences in morphology and depth from the lower to the 
upper estuary, play a key role in setting the circulation patterns 
found in the system. This influence is observed in major and 
shallow water tidal constituents, tidal asymmetry, tidal dissipation 
and tidal ellipses patterns. 

Results of amplitude and form factor have given evidence that 
tide is semidiurnal. The M2 amphtude increases within the estuary 
from Cascais to Ponta da Erva. Even though losing energy to the 
shallow-water constituents, the semidiurnal tide grows up as it 
propagates from the mouth until Ponta da Erva, decreasing 
further upstream at Vila Franca de Xira. This pattern suggests that 
between Ponta da Erva and Vila Franca de Xira there is a fraction 
of the tidal energy reflected somewhere, which is consistent with 
the existence of a resonant mode with the Afj period. 

The results also showed differences in the tide propagation in 
two main areas of the estuary (the upper and lower estuary), 
resulting from a balance between convergence/divergence and 
bottom friction and advective effects. There is a phase delay as the 
tidal wave moves toward the inshore zone of the system. In the 
shallow areas the is strongly attenuated, while the opposite 
occurs for the M4 and Mg, which are strongly amplified, thus 
inducing tidal asymmetry. The amplitudes of both shallow-water 
constituents M4 and Mg enlarge along the estuary, from the mouth 
to the head. The M4 constituent amplitude increases from Cascais 
to Cabo Ruivo, whereas Mg amplitude is smoothly enlarged from 
the mouth to Cabo Ruivo, and with a higher gradient until de 
estuary head. Considering that the M4 generation is essentially due 
to the nonlinear effects of advection and that Mg is mainly induced 
by bottom friction, these patt(-rns suggest that the advective 
accelerations are dominant from tiic estuary mouth to Cabo 
Ruivo, while bottom friction takes a larger significance upstream 
that zone. 

The ellipses orientation reveals the importance of the bottom 
topography and of the bay geometry on the tidal fluxes, with the 
main flux through the deeper zones. Tidal asymmetry parameters 
indicate ebb dominance for most of the estuary, so there is a general 
trend for a sediment exportation, and consequentiy for erosion 
with the natural deepening of these areas. The results for the tidal 
dissipation showed that the dissipation is higher in deeper areas 
where tidal currents are stronger, which is the case of the main 
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channel (upper estuary) and in two preferential channels of the 
bay. The analysis of this parameter simultaneously with tidal 
ellipses orientation show that the flow in the Tagus estuary occurs 
by preferential channels, where current speeds are higher due to 
their highest depths. 
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